#/usr/bin/perl                                                                                                                           

# GLOBAL

my $grid_size_x = 300;
my $grid_size_y = 300;
my $locitot = 99;
my $name = "sim";

# Seed the random number generator
srand (0);

sub draw_rand{
    my ($n,$tot) = @_;
    my %h, @arr;
    while(keys(%h)<$n){
        $r = int(rand($tot));
        if(!$h{$r}){
            $h{$r}=1;
            push(@arr, $r);
        }
    }
    @arr = sort {$a <=> $b} @arr;
    return @arr;
}

# draw_rand_hash(n,tot) draws n random numbers from 0..(tot-1) without duplicates 
sub draw_rand_hash{
    my ($n,$tot) = @_;
    my %h;
    while(keys(%h)<$n){
        $r = int(rand($tot));
        if(!$h{$r}){
            $h{$r}=1;
        }
    }
    return %h;
}

# Draw random x cells, but since we have hex grid
# x = 0, 0.5, 1, 1.5 ... 511, 511.5
# No two cells have the same coordinates
sub draw_rand_cells{
    my ($n,$tot) = @_;
    my %h, $size=0;
    while($size<$n){
        $x = int(rand($tot));
	# either add 0.5 to both, or not
	$offset = 0.5*int(rand(2));
	$x+=$offset;
	if($offset<0.2){
	    $y = int(rand($tot/2))*1.5;	    
	} else {
	    $y = int(rand($tot/2))*1.5+0.75;	    
	}
        if(!$h{$x}{$y}){
            $h{$x}{$y}=1;
	    $size++;
        }
    }
    return %h;
}

# Start from sweep 1
# do not use my, so that these hashes are global
 %grid;
 %patterns;
 %indexx, %indexy; 
for $sw ($ARGV[0]..$ARGV[1]){

    my %grid_coords_x =();
    my %grid_coords_y =();
    open (in,"$name$sw/$name$sw.crypts");
    while(<in>){
	chomp;
	@row=split(/\,/);
	$grid_coords_x{$row[0]}=$row[1];
	$grid_coords_y{$row[0]}=$row[2];
    }
    close in;

    %patterns=();
    # Important = dead crypts have cloneId=-1, so include cloneId=-1 in patterns
    # the pattern of a dead crypt is wild type...
    # Try hashes of arrays to prevent memory explosions
    #my @temp;
    #for my $loc (0..243){
#	$temp[$loc]=100;
 #   }
    # load up an array
    #$patterns{-1}=[@temp];
    # Do not include dead cells in the divergence computation
    # is the correct way to do it
    #$patterns{-1}=1;

    # Load all patterns in memory
    open (in,"$name$sw/$name$sw.patterns");
    while(<in>){
	chomp;
	@row=split(/\,/);
	$patterns{$row[0]}="";
	for my $loc (105..109){
	    # Hash only the locations that are different from 
	    # the wild type "0" microsattelite length	    
	    # Shaka Zulu proposes hashing only the clone_id key
	    # and as value concatenate "locus,value,locus,value"
	    if($row[$loc+1]!=0){
		# Leading comma ",loc,val,loc,val"
		$patterns{$row[0]}=$patterns{$row[0]}.",".$loc.",".$row[$loc+1];
	    }
	}
    }
    close in;

    print STDERR "read done $sw\n";

    my $day = 182; 
    my $generation;
    my $firstline=1;
    # Go through the log file of grid state
    %grid=();
    open (i1,"$name$sw/$name$sw.log");
    open (out, ">$name$sw/$name$sw.hit.log");
    while(<i1>){
	chomp;
	@row=split(/\,/);
	if($firstline){
            $generation=$row[3];
            $firstline=0;
        }
	if($row[2]!=$day){
	    
	    # compute divergence
	    #print STDERR "start comp diverg\n";
	    # Compute global divergence for $n cells
	    my $n = 300;
	    
	    # Draw cells
	    %hashxy = ();
	    %hashxy = draw_rand_cells($n,$grid_size_x);
	    %indexx = ();
	    %indexy = ();
	    my $i=0;
	    for my $x (keys %hashxy) {
		for my $y (keys %{$hashxy{$x}}) {
		    $indexx{$i}=$x;
		    $indexy{$i}=$y;
		    $i++;
		}
	    }
	    
	    my @count = compute_percentages($n,\%grid,\%patterns,\%indexx,\%indexy);
	    print out "$day,$generation,";
	    for my $p (0..4){
		if($p<4){
		    print out sprintf("%.3f", $count[$p]/$n).",";
		} else {
		    print out sprintf("%.3f", $count[$p]/$n)."\n";
		}		    
	    }

	    # increase day
	    if($day==7280){
		$day=7300;
	    } else {
		$day+=182;
	    }
	    $generation=$row[3];
	    %grid=();
	    # Important - throw it as a double in hash
	    $grid{$grid_coords_x{$row[0]}*1.0}{$grid_coords_y{$row[0]}*1.0}=$row[1];
	} else {
	    # grid X coordinate, Y coordinate, clone_id
	    $grid{$grid_coords_x{$row[0]}*1.0}{$grid_coords_y{$row[0]}*1.0}=$row[1];
	    #print STDERR $grid_coords_x{$row[0]}." ".$grid_coords_y{$row[0]}." ".$row[1]."\n";
	}
    }
    close i1;    

    #print STDERR "loop done\n";

    # For day 7300
    # Compute global divergence for $n cells
    my $n = 300;
    
    # Draw cells
    %hashxy = ();
    %hashxy = draw_rand_cells($n,$grid_size_x);
    %indexx = ();
    %indexy = ();
    my $i=0;
    for my $x (keys %hashxy) {
	for my $y (keys %{$hashxy{$x}}) {
	    $indexx{$i}=$x;
	    $indexy{$i}=$y;
	    $i++;
	}
    }
    
    my @count = compute_percentages($n,\%grid,\%patterns,\%indexx,\%indexy);
    print out "$day,$generation,";
    for my $p (0..4){
	if($p<4){
	    print out sprintf("%.3f", $count[$p]/$n).",";
	} else {
	    print out sprintf("%.3f", $count[$p]/$n)."\n";
	}		    
    }
    
    close out;    
    print STDERR "computing perc done $sw\n";
}



sub compute_percentages{
    # Number of cells
    my ($n,$gridRef,$patternsRef,$indexxRef,$indexyRef) = @_;
    my %grid = %$gridRef;
    my %patterns = %$patternsRef;
    my %indexx = %$indexxRef;
    my %indexy = %$indexyRef;

    my @count;
    for my $i (0..4){
	$count[$i] = 0;
    }
    my $dead_cells_count = 0;
    # Compare pattern of cell i to patterns of cell j
    # A total of n*(n-1)/2 comparisons 
    for my $i (0..($n-1)){
	my @col;
	# get the clone_id
	my $cid = $grid{$indexx{$i}}{$indexy{$i}};
	if(length($patterns{$cid})>0){
	    my @cell_one = split(/\,/,$patterns{$cid});
	    my $hit_count=$#cell_one/2;
	    $count[$hit_count-1]++;	
	}
    }
    
    return @count;
}
